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The effects of opacity on gravitational stability in protoplanetary 
discs 
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ABSTRACT 

In this paper we consider the effects of opacity regimes on the stability of self-gravitating 
protoplanetary discs to fragmentation into bound objects. Using a self-consistent 1-D viscous 
disc model, we show that the ratio of local cooling to dynamical timescales Q? coo i has a strong 
dependence on the local temperature. We investigate the effects of temperature-dependent 
cooling functions on the disc gravitational stability through controlled numerical experiments 
using an SPH code. We find that such cooling functions raise the susceptibility of discs to frag- 
mentation through the influence of temperature perturbations - the average value of Qf coo i has 
to increase to prevent local variability leading to collapse. We find the effects of temperature- 
dependence to be most significant in the 'opacity gap' associated with dust sublimation, where 
the average value of Qf coo i at fragmentation is increased by over an order of magnitude. We 
then use this result to predict where protoplanetary discs will fragment into bound objects, 
in terms of radius and accretion rate. We find that without temperature dependence, for radii 
< 10 AU a very large accretion rate ~ 10" 3 M Q yr" 1 is required for fragmentation, but that 
this is reduced to 10~ 4 M G yr -1 with temperature-dependent cooling. We also find that the 
stability of discs with accretion rates < 10~ 7 M G yr~' at radii > 50 AU is enhanced by a lower 
background temperature if the disc becomes optically thin. 
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1 INTRODUCTION 

The formation of planets within protoplanetary discs is a sub- 
ject that attracts considerable interest, with two main compet- 
ing schools of thought. The core accretion-gas capt ure model 
dLissauer. 19 93 Lissauer & Stevenson. 200jKlahr. 2 008 ) posits hi- 
erarchical growth, with the collisional coagulation of dust grains 
initially leading to centimetre-sized particles, and thence on to 
planetesimals and rocky planets. Once a critical mass is reached, 
it is then possible to accrete a gaseous envelope and hence 
form giant Jupiter-like planets. Various observations have suc- 
cessfully confirme d this mode of planet formation, for exam ple 
iMarcv et al. (2005I) . lDodson-Robinson & Bodenheimer (2009t) . 

However, this model cannot explain all the available observa- 
tions . iKennedv & Kenvon (2008h show that beyond approximately 
20 AU, the timescales for giant planet formation via core accretion 
exceed the expected disc lifetime of approximately 10 Myr, imply- 
ing that no planets should be detected in this region. However, re- 
cent observations of HR8799 with the Keck and Gemini telescopes 
have produced direct images of giant planets (5 - 13 Mj) orbiting at 



radii of up to ~ 70 AU dMarois et al.. 2008|). Similar observations 
of ot her systems (e.g. / ? Pic b dLagrange et al.. 2009h and Formal- 
huat dKalas et al.. 2008|)) and theor etical work on the formation of 
2MASS 1207b dLodato et al„ 20051) have suggested that there is an- 
other mechanism for planet formation at work, and this is thought 
to be the effect of gravitational instabilities within the protoplane- 
tary discs themselves. 

In protoplanetary discs where the self-gravity of the 
gas is dynamically important, direct gravitational collapse 
of locally Jeans-unstable over-densiti es within the disc 
dBoss. 199lBoss. 1998fc>urisen et al.. 20071) would also pro- 
duce giant planets very rapidly, on the local dynamical 
timescale. A similar process of gravitational instability lead- 
ing to local collapse is a strong candidate for the forma- 
tion of stellar discs around Active Galactic Nuclei (AGN) 
dNavakshin et al.. 2007) and those observed in our own Galacti c 



Centre dLevin & Beloborodov. 2003lNavakshin & Cuadra! 20051) . 
and in the context of protostellar discs may furthermore be 
responsible for the for mation of brown dwarves and other low 
mass stellar companions dStamatellos et al.. 2007aT) . 
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The emergence of the gravitational instability within a disc is 
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governed by the parameter O jToomre. 1964) . which for a gaseous 
Keplerian disc is given by 



(i) 



This encapsulates the balance between the stabilising effects of 
rotation {£l(R) is the angular frequency at radius R) and ther- 
mal pressure (c B (R) is the sound speed) and the destabilising ef- 
fect of the disc self-gravity via the surface density l.(R). When 
Q < 1 the instability is initiated, leading to the presence of spi- 
ral density waves within the disc which, depending on the lo- 
cal cooling rate, may persist in a sel f-regulated quasi-stable state 
jGammie~ 00 1 Lodato & Rice, 2004) or may fragm ent into bound 
clumps jjohnson & Gammie. 200 jRice et al.. 20051) . Discs that are 
sufficiently cool are therefore expected to be susceptible to the 
gravitational instability, and it is thought that at least in the early 
stages of stellar evolut ion, many discs enter the self-gravitating 
phase l lHartmann. 20091) . 

Once the gravitational instability is initiated, heat is input to 
the disc on the dynamical timescale through th e passage of spiral 
compression/shock waves l lCossins et al.. 20091) . Various numerical 
studies using both 2D and 3D models of self-gravitating discs have 
produced the result that, in order to induce fragmentation, the disc 
must be able to cool on a ti mescale faster than a few times th e local 
dynamical time, /d yn = fl 1 jGammie, 200llRice et al. 20051) . This 
condition is likely to occur only at relatively large radii (~ 100 AU), 
on the ass umption that stellar or external irradiation of the disc is 
negligible dRafikov. 2009lStamatellos & Whitworth. 2009bh . 

These models have generally used a cooling rate prescribed by 
using a fixed ratio between the local cooling (/ coo i) and dynamical 
(Or 1 ) times, such that 



(2) 



for some constant ft throughout the radial extent of the disc. 
Various authors (e.g. iGammie 200llRice et al. 20051) have found 
that fragmentation occurs whenever flf CO oi ~ 3 - 7. By using 
a more realistic cooling f ramework based on the optical depth, 
Ijohnson & Gammie (2003b found that the fragmentation bound- 
ary (defined hereafter as the ratio of the cooling to dynamical 
timescales, flf coo i at fragmentation) may in fact be over an order 
of magnitude greater than this, leading to an enhanced tendency 
towards fragmentation. This variation in £lf coo i they ascribed to the 
implicit dependence of the cooling function on the disc opacity, and 
hence on temperature. 

Using the opacity tables of iBell & Lin (1994) it is clear that 
the opacity is a strong function of temperature in certain regimes, 
and by modelling protoplanetary discs as optically thick in the 
Rosseland mean sense, flf C ooi shows power law dependencies on 
both the local temperature and density. In cases where this depen- 
dence is strong, it is therefore possible that small temperature fluc- 
tuations may push the local value of fl/ CO oi below the fragmentation 
boundary, even when the average value is significantly above it. 

In this paper we therefore seek to investigate and clarify 
the exact relationship between the fragmentation boundary 
and the temperature dependence of flf CO oi, using a Smoothed 
Particle Hydrodynamics (sph) code to conduct global, 3D nu- 
merical simulations of discs where the cooling time follows 
a power-law dependence on the local temperature. In addi- 
tion, various studies have shown that in a quasi-steady state 
the gravitational instability may be modelled pseudo-viscously 
dLodato & Rice 2005fcossins et al. 200 ^Clarke 200jRafikov 20091) . 



We therefore use the a-prescription of Shakura & Sunvaev (19731) 
and the assumption of local thermal equilibrium, where 



fifco 



1 



9 7(7 - l)ff 



(3) 



and 7 is the ratio of specific heats, to construct an analytical model 
of the opacity regimes present within a marginally gravitationally 
stable disc. From this we can therefore predict analytically if and 
where such discs would become prone to fragmentation, and also 
compare these results to those from the more co mplex simula- 
tio ns where radiative transfer is mod elled, such as iBolev (20091) 
and lStamatellos & Whitworth (2009bh . 

The structure of this paper is therefore as follows. In Section 2 
we discuss some of the theoretical results relevant to protoplanetary 
discs, and introduce a simplified cooling function derived from the 
various opacity regimes. We further consider the effects we expect 
these cooling prescriptions to have on the susceptibility of proto- 
planetary discs to fragmentation. In Section 3 we briefly outline the 
numerical modelling techniques used in our simulations and detail 
our initial conditions. In Section 4 we present the results from these 
simulations, before proceeding to collate these with the analytical 
predictions in Section 5. Finally in Section 6 we discuss the rami- 
fications of our work and the conclusions that may be drawn from 
it. 



2 THEORETICAL RESULTS 

In this section we derive analytical results for the dependence of 
the cooling timescale / CO oi on temperature and density, such as we 
might expect to find in a quasi-gravitationally stable protoplanetary 
disc environment. We also consider analytically the effects that a 
(specifically) temperature dependent cooling time will have on the 
stability of such a disc to fragmentation. 



2.1 flf coo i in the Optically Thick Regime 

As in the case of I Gammie" (200 ll) . we may start from the following 
basic equations: 

t/Z 

fcool = -T-, (4) 

A 

T * pH K , (5) 
Z = 2pH, (6) 
YRT 

C; = — 

where U is the specific internal energy, A is the cooling rate per unit 
area, r is the optical depth, p is the (volume) density, H = c s /fl is 
the disc scale height, k is the opacity, 7 is the ratio of specific heats, 
ft = k/mn is the universal gas constant (k being the Boltzmann 
constant and mu the mass of a hydrogen atom), T is the local mid- 
plane temperature and yu is the mean molecular weight of the gas. 
Note that the factor of two in equation [6] arises from there being 
two faces of the disc from which to radiate. 

In the case where the disc is optically thick (in terms of the 
Rosseland mean), then the cooling rate per surface area A may be 
given as 



(7) 



A : 



I60T 

3t 



4 



(8) 



where cr is the Stefan-Boltzmann constant. We note that this is 
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strictly valid only in the case where energy is transported radia- 
tively within the disc — convective transport or stra tification withi n 
the disc will alter this relationship (see for example iRafikov 20071) . 
For the purely radiative case, the vertical temperature structure 
of the disc is therefore accounted for via this formalism, and 
is characterised by the midplane temperature T and the optical 
depth t. In order to prevent divergence of this cooling function 
at low optical depths and to inte rpolate smoothly into the opti- 
cally thin regime, others including I Johnson & Gammie (20031) and 
iRice & Armitage (20091) have used a cooling function of the form 



16o-r 4 



(9) 



which becomes directly proportional to the optical depth in the op- 
tically thin limit. In general however, we find that discs only be- 
come optically thin at large radii, and that this correction is there- 
fore only relevant to the case where the cooling is dominated by 
ices. 

Furthermore, we note that for systems where the stellar mass 
dominates over that of the disc, the density p may be approximated 
by 

M, 

9 ~ 2ttR 3 Q 



(10) 



where M, is the mass of the central star and R the radial distance 
from the central star, and therefore we have Q. 2 = 2nGpQ in the 
case of Keplerian rotation, with G being the universal gravitation 
constant. Recalling also that c 2 = Uy(y - 1), equations I4l-I10lmay 
be rearranged to show that in the optically thick case, the ratio of 
cooling to dynamical times should be 



3<R 2 



y 



Q- i,2 P 3 > 2 T- 2 . 



(11) 



8cr V2?rG 7 ~ 1 P 2 

iBell & Lin (19941) found that the opacity can be reasonably 
well approximated by power-law dependencies on temperature and 
density, such that 



K p 



a jb 



(12) 



Specific values of a, b and kq apply for each opacity regime, such 
that the value of k varies continuously over the regime boundaries. 
Using these approximations, we find that the Q/ coo i value for the 
various opacity regimes can be given by 



Clt a 



3<R 2 



7*0 



1) 



q-I/2 ^"+3/2 T b 



(13) 



For each opacity regime, the constant k , the exponents a and b, 
the transition temperatures between the regimes and the functional 
dependence of Q.t coo \ on temperature and density are given in Ta- 
ble [T] It should be noted that for the purposes of these tables the 
temperature and density should be measured in cgs units. 



2.2 Effects of Temperature Dependence on Fragmentation 

We now specifically consider the effects of temperature fluctuations 
on the stability of a disc to fragmentation, using a simplified cool- 
ing prescription derived from a consideration of equation [T3l 

In the previous section it was noted that the ratio of the lo- 
cal cooling and dynamical times f2/ c00 i has a direct dependence 
on the local mid-plane temperature T. Given that (from Table Q} 
this dependence is generally much stronger than that on density, it 
is physically reasonable to consider a simplified cooling function 



where we only include the effects of temperature, and where we 
define the cooling time via the relationship 



■'(?)" 



(14) 



for some general value of the cooling exponent n and cooling pa- 
rameter p. Here T is the azimuthally averaged mid-plane tempera- 
ture T in thermal equilibrium, and thus we see that when thermal 
equilibrium is reached, the average cooling timescale is expected 
to reduce to (f2f cool ) a /?, with a fragmentation boundary /?„ as- 
sociated with each value of n. In particular, wit h n = 0, at frag - 
mentation we ha ve Q.t coo i = (Q.t coo i) = f} , which iGammie (200 lh . 
IRice et al. (20051) and others have found to be in the range 3-7. 

In the case of temperature dependent cooling (where n t 0), 
if the equilibrium value of the cooling parameter p > f3 , the disc 
may still fragment due to temperature fluctuations leading to a short 
term (relative to the dynamical timescale) decrease in the instan- 
taneous value of {S to less than the threshold value. For a power- 
law index n, in order to calculate the value /?„ of the equilibrium 
cooling parameter below which fragmentation occurs, we make the 
assumption that fragmentation takes place wherever the instanta- 
neous value of f2?cooi is held at or below the critical value /Jo for 
longer than a dynamical time, independent of the mechanism by 
which the cooling is effected. If we therefore consider temperature 
fluctuations such that T = f + 6T, we find that at the fragmentation 
boundary 



■*(«♦?) 



(15) 



In lCossins et al. (20091) for the case where M iix /M, = 0.1 we 
found that on average the strength of the surface density perturba- 
tions 61, /Y, can be linked to the strength of the cooling through the 
following relationship. 



2 J ~ (Qt coal ) [/2 ' 



(16) 



where angle brackets denote the RMS value. In a similar manner 
we may say that 



(-) = 



(17) 



where k is to be defined empirically. At fragmentation therefore we 
have 



(f) 



k 



(18) 



noting that by construction for a given index n, at fragmentation 
(f2f coo i> = /?„. Combining this with equation Q3] we find that in the 
case where the cooling is allowed to vary with temperature as per 
equation [14] the fragmentation boundary /?„ satisfies the following 
equation; 



Pn 



(19) 



This implicit equation can therefore be solved to find the value of 
the fragmentation boundary fS„ for all n > -2 (below this /3„ be- 
comes undefined), as shown later in Table|4] 
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Table 1. Details of the various optical regimes by type, showing the transition temp eratures and the fu nctional dependence of Clt coo \ on the temperature and 
density in the optically thick regime. Note that all values are quoted in cgs units. See lBell & Lin (19941) for further details. 



„ .. „ . , 9 _i. . Temperature Range (K) _ , 

Opacity Regime *:o (cm- g ) a b t Dependence ot 12^ 



Ices 


2 X 1(T 4 





2 





166.810 


p 3/2 


Sublimation of Ices 


2x 10 16 





-7 


166.810 


202.677 


p 3/2 T -9 


Dust Grains 


1 X 10-' 





1/2 


202.677 


2286.77 p 2/49 


p 3/2 j-5/2 


Sublimation of Dust Grains 


2x 10 81 


1 


-24 


2286.77 p 21 * 9 


2029.76 p 1/81 


p 5/2 j-26 


Molecules 


1 X 1(T 8 


2/3 


3 


2029.76 p 1/81 


10000.0 p 1/21 


p 13/6 jl 


Hydrogen scattering 


1 x 10~ 36 


1/3 


10 


10000.0 p 1/21 


31 195.2 p 4/75 


pi 1/6 ji 


Bound-Free & Free-Free 


1.5 x 10 20 


1 


-5/2 


31 195.2 p 4/75 


1.79393 x 10 8 p 2 ' 5 


p S/2 T -9/2 


Electron scattering 


0.348 








1.79393 x 10 8 p 2/5 




p 3/2 j-2 



3 NUMERICAL SET UP 
3.1 The SPH code 

All of the simulations presented hereafter were performed us- 
ing a 3D smoothed particle hydrodynamics (sph) code, a La- 
grangian hydrody namics code capable of mode lling self-gravity 
(see for example, iBenz 1990l iMonaghan 1992h . The code self- 
consistently incorporates the so-ca lled V/i terms to ensure en - 
ergy conservation, as de scribed in ISpringel & Hernquist (20021) . 
IPrice & Monaghan (20071) . All particles evolve according to indi- 
vidual time-steps govern ed by the Courant c ondition, a force con- 
dition (Monagha n. 19921) . an integrator limit dBate et al.. 19951) and 
an additional condition that ensures the local timestep is always less 
than the local cooling time. 

We have modelled our systems as a single point mass 
(onto which gas particles may accrete if they enter within a 
giv en sink radius a nd satisfy certain boundness conditions — 
see iBate et al. 19951) orbited by 500,000 sph gas particles; a set 
up common to many other sph simulations of such systems (e.g. 
iRice et al. (200jLodato & Rice (2004Lodato & Rice (200 jClarke et al 
but at a higher resolution than most. The central object is free to 
move under the gravitational influence of the disc. 

In common with many other simula- 
tions where cooling is being investigated 

dGammie (200llLodato & Rice (2005fcossins et al. (20091) for 
example) we use a simple implementation of the following form; 

dui u: 

~T = ( 20 ) 

«f 'cool,; 

where z<, and f cool / are the specific internal energy and cooling time 
associated with each particle respectively. The cooling time is al- 
lowed to vary with the particle temperature T t in such a manner 
that 

CVcooi, =/3(j) . (2D 

where f2, is the angular velocity of the particle, f is the equilibrium 
temperature, and fi and n are input values held constant throughout 
any given simulation. Given that T ~ c 2 , equation [T] shows that for 
a given value of the surface density S this is equivalent to 

JVcooy=$(|fj , (22) 

where again Q, is the value of the Q parameter evaluated at each 
particle, and Q is the expected equilibrium value of Q, which we 
take to be 1 throughout. Note that a priori we do not know exactly 
what the equilibrium value of Q will be once the gravitational in- 



stability has saturated. Indeed as we shall see this turns out to be 
slightly greater than unity, but still such that Q « 1. The effective 
value of the cooling parameter is given by 

/3=PQ- 2 ", (23) 

where Q is the actual value to which the simulations settle. Since 
we are exploring relatively large values of n, p can vary signifi- 
cantly from our input value for even small changes in Q. 

Finally we calculate the equivalent surface density E, (and thus 
Qi) at the radial location of each particle R t by dividing up the disc 
into (cylindrical) annuli, calculating the surface density for each 
annulus, and then interpolating radially to obtain To pre- 

vent boundary effects, for simulations where n > 1.0 the tempera- 
ture dependent effects are limited to an annulus 15 < R < 20 (in 
code units — note that initially R m = 0.25 and i? out = 25.0). At 
other radii we keep f2/ coo i = 8, a value chosen to suppress fragmen- 
tation in regions outs ide the annulus of interest (see for instance 
lAlexander et al. 20081) . 

All simulations hav e been run with the particles modelled as 
(2fla?ff@asiga«tvdtlG2<MJ^itio of specific heats y = 5/3. Heat ad- 
dition is allowed for via PdV work and shock heating. Artificial 
viscosity has been included through the standard sph formalism, 
with ffspH = 0.1 and /3 SPH = 0.2 — although these values are 
smaller than those commonly used in sph simulations, this limits 
th e transport and heating induced by artificial viscosity. As shown 
in lLodato & Rice (20041) . with this choice of parameters the trans- 
port of energy and angular momentum due to artificial viscosity is 
a factor of 10 smaller than that due to gravitational perturbations, 
while we are still able to resolve the weak shocks occurring in our 
simulations. 

By using the cooling prescription outlined above in equa- 
tion [22] the rate at which the disc cools is governed by the dimen- 
sionless parameters Q, ft and n, and the cooling is thereby imple- 
mented scale free. The governing equations of the entire simulation 
can therefore likewise be recast in dimensionless form. In common 
with the previous sph simulations mentioned above, we define the 
unit mass to be that of the central object - the total disc mass and 
individual particle masses are therefore expressed as fractions of 
the central object mass. We can self-consistently define an arbitrary 
unit (cylindrical) radius Rq, and thus, with G = 1, the unit timestep 
is the dynamical time f dyn = Q. [ at radius R = 1. 

3.2 Initial conditions 

All our simulations model a central object of mass M„, surrounded 
by a gaseous disc of mass Mdi sc = 0.1 M,. We have used an ini- 
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Table 2. Table of simulations run for various values of the cooling exponent 
n and rate fi. Note that since many of these simulations were run concur- 
rently, there is a degree of overlap in the /J values used. 



Exponent (n) 



Input cooling parameter 0) 



0.0 3, 4, 4.5. 5, 6 

0.5 4,4.5,5,5.5,6 

1.0 3,4,5,6,7,8,9,10 

1.5 7,8,9,10,11 

2.0 10,11,12,13,14,15,16,17,18 

3.0 20, 22.5, 25, 27.5, 30, 32.5, 35, 37.5, 40 



tial surface density profile S ~ R~ 3/2 , which implies that in the 
marginally stable state where Q » 1, the disc temperature profile 
should be approximately flat for a Keplerian rotation curve. Since 
the surface density evolves on the viscous time f v i sc » fdyn = ^ 1 
this profile remains roughly unchanged throughout the simulations. 
Radially the disc extends from i?j n = 0.25 to R out = 25.0, as mea- 
sured in the code units described above. The disc is initially in 
approximate hydrostatic equilibrium in a Gaussian distribution of 
particles with scale height H. The azimuthal velo cities take into ac- 
count both a pressure correction l lLodato. 20071) and the enclosed 
disc mass. In both cases, any variation from dynamical equilibrium 
is washed out on the dynamical timescale. 

The initial temperature profile is c 2 ~ R~ [/2 and is such that 
the minimum value of the Toomre parameter <2 rain = 2 occurs at 
the outer edge of the disc. In this manner the disc is initially grav- 
itationally stable throughout. Note that the disc is not initially in 
thermal equilibrium - heat is not input to the disc until gravitational 
instabilities are initiated. 



3.3 Simulations run 

Since our simulations use a slightly different surface density pro- 
file to that used by previous authors (L ~ R ~ 3/2 , cf. S ~ R' 1 in 
iRice et al. 2001 S ~ R ~ 7/4 in lRice et al. 20031) we initially ran five 
simulations at various values of p with the cooling exponent n set 
equal to zero, to find the fragmentation boundary in the case where 
the cooling is independent of temperature. Thereafter, simulations 
were run at various /? values as n was incremented up to n = 3, to 
ascertain the fragmentation boundary in each case. A summary of 
all the simulations run is given in Table|2] 



4 SIMULATION RESULTS 
4.1 Detecting Fragmentation 

First of all it is useful to explain how fragmentation has been de- 
tected in our simulations. Throughout all the numerical simula- 
tions run, the maximum density over all particles has been tracked 
as a function of elapsed time. In the case of a non-fragmenting 
disc, the maximum always occurs at the inner edge of the disc 
(as would be expected), and is relatively stable over time. How- 
ever, once a fragment forms, this maximum density (now corre- 
sponding to the radius at which the fragment forms) rises expo- 
nentially, on its own dynamical timescale. An example is shown in 
Fig. U] and the various changes in gradient correspond to various 
fragments at different radii (and thus with differing growth rates) 
achieving peak density. A simil ar increase in the central density o f 
proto-fragments is observed in lStamatellos & Whitworth (2009ah . 
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Figure 1. Maximum density plot showing the characteristic rise due to frag- 
ment formation, seen here for the simulation where ft = 4.0, n = (where 
the cooling is independent of temperature). There is clear evidence of frag- 
ment formation at t as 20, 000, with both density and time being shown in 
code units. 



although the timescales differ due to the use of different equations 
of state. 

This rise in the maximum density has therefore been used 
throughout as a tracer of fragment formation, and the evolution has 
been followed until the fragments are at least four orders of magni- 
tude greater than the original peak density. 



4.2 Averaging Techniques 

Throughout the following analysis, we have defined the average 
value of a (strictly positive) quantity, which we denote by an over- 
bar, as the geometric mean of the particle quantities. The reason 
for this is we find that in the "gravo-turbulent" equilibrium state 
properties such as the temperature, density and Q value are log- 
normally distributed. This is shown for example in Fig.[2] where the 
temperature data from the simulation match a predicted log-normal 
distribution to within one percent. (Note the reduced radial range 
to reduce the effect of the inherent gradual reduction in tempera- 
ture with radius.) The geometric mean being precisely equivalent 
to the exponential of the arithmetic mean of the logged values, this 
process recovers the mean value of the normal distribution of In T. 

Similarly, to calculate the perturbation strengths (e.g 6A/A for 
some quantity A) we note that 



SA 

T 



dA 
~A 



din A. 



(24) 



The RMS value of SA/A is then equivalent to the standard deviation 
of In A, which again can be recovered directly from the log-normal 
distribution. Referring again to Fig. [2] we therefore see that T = 
10 -4.498 _ 3 1?7 x 1Q -5 ^ code units)j and that ST /f = (T = o.348. 
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LogN(ln T, a 2 ) = Log_\(-4.498, 0.121) 
Simulation / \" 
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Figure 2. Distribution of particle temperatures for 16.25 < R < 18.75, and 
a predicted log-normal distribution based on the same data. The two are 
equal to within approximately 1%. 



Figure 3. Plot of Q against radius for various values of in the temperature 
independent case n = 0. For the fragmenting cases (8 < 4.25) the values 
shown are from immediately prior to fragmentation. 



4.3 Equilibrium States 

First of all, we need to determine the exact value of the fragmen- 
tation boundary in the case where n = (and thus where p = /?), 
which we denote by /Jo- As seen in Table [2] simulations were run 
at various values 3.0 < f3 < 6.0, and we find that the boundary lies 
between 4.0 and 4.5. We therefore take the critical value as being 
the midpoint, such that {S = 4.25. 

Continuing with the n = case, we find throughout that the 
value of Q to which the simulations settle is slightly above unity. 
The steady state values (time averaged over 1000 timesteps) are 
shown for various ft values in Fig. [3] and we see that the average Q 
value is approximately 1.091, where we have averaged over both /? 
and radius (where 15 < R < 20, for comparison with simulations 
with higher n). We further note that there is scatter of ~ 10% about 
this average, and (although not shown) this is equally true of the 
simulations where n > 0. 

Note then that for large n the effective value of the cooling pa- 
rameter p at any given radius may be substantially different from 
the numerical input value ft = f3Q ln (see equation [23} that we use 
to characterise the cooling law. In order to determine the fragmen- 
tation boundary with any accuracy, we therefore need to consider 
the true value of f3 rather than the input value 0. 

4.4 Cooling Strength and Temperature Fluctuations 

In order to characterise the fragmentation boundary, it is necessary 
that we validate the assumption encompassed by equation [17] that 
the temperature perturbation strength is correlated to that of the 
applied cooling. Using the method outlined above in section |4~2l 
for each simulation we can calculate azimuthally averaged RMS 
values for the strength of the temperature fluctuations, which we 
denote by (ST /T ^. Where n = 0, these temperature perturbations 
are plotted as a function of radius for various values of p in Fig. [4] 
where we see that there is a systematic decrease in the perturbation 



0.8 , , , , | , , , , | , , , , | , , , , | , , , , | , , , 




R [R„] 

Figure 4. Plot showing the strength of temperature perturbations within the 
disc as a function of radius and fi for the temperature independent case, 
where n = 0. 

strength with increasing /?, and also that the perturbation strength is 
almost constant with radius across the self-regulating region (5 < 
R < 25) of the disc. Using equation [T7] we can therefore calculate 
an empirical value for k, and hence averaging both radially (for 
15 < R < 20 as before) and over the available values of ft we find 
k = 1.170 where n = 0. 

Furthermore, we note that in the temperature dependent case 
(where n + 0), by construction the average value (f2f cool ) is simply 
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Table 3. Table showing the fragmentation boundaries obtained from the 
simulations. The central columns show respectively the highest fragment- 
ing and lowest non-fragmenting values of /? simulated, with /}„ being the 
midpoint of these. Throughout, /J is calculated using equation l23l 



Exponent (n) 


Effective cooling rate (Ji) 
Fragmenting Non-Fragmenting 


A, 


0.0 


4.000 


4.500 


4.250 


0.5 


4.825 


5.263 


5.044 


1.0 


5.915 


6.654 


6.284 


1.5 


6.949 


7.644 


7.296 


2.0 


8.458 


9.022 


8.740 


3.0 


10.051 


11.056 


10.554 



the effective value of the cooling strength, /3. We can therefore cal- 
culate the value of k for cases where n t 0, and we find that again k 
remains constant both with the index n and with radius. Hence we 
take the value of k to be 1.170, as in the n = case, and empirically 
we may therefore say that on average 

for all n. 



(25) 



4.5 The Fragmentation Boundary 

We are now in a position to predict empirically the fragmentation 
boundary in the case where n t 0, and to compare this directly 
with the results of our simulations. Table[3]shows the fragmentation 
boundary /?„ as obtained from our simulations, where once again it 
is taken as the average of the highest fragmenting and lowest non- 
fragmenting values of /? simulated. We find that as expected, there 
is indeed a rise in the fragmentation boundary as the dependence of 
the cooling on temperature increases. This variation of the fragmen- 
tation boundary is shown against the cooling exponent n in Fig. [5] 
(where the error bars show the upper and lower bounds from Ta- 
ble [3} along with predicted values generated using the following 
empirically defined implicit relationship 



1 + 



1.170 



(26) 



where we have used /?o = 4.25. Clear from this plot is the fact that 
the predictions are a very good match to the data observed, and our 
theoretical model, in which the increased tendency for fragmenta- 
tion is due to the effects of temperature fluctuations on the cooling 
rate, is therefore valid. The transition zone shown is bounded by 
curves corresponding to predictions using /? = 4.00 and 4.50, the 
upper and lower bounds for /?o we obtained through our simula- 
tions. 

4.6 Statistical Analysis 

The effects of temperature perturbations on the fragmentation 
boundary can be neatly illustrated statistically, if we assume that 
the distribution of temperatures about the geometric mean In f is 
log-Normal (as found in our simulations). Using standard notation 
we can therefore say that 



In T ~ N(\n T, a 1 ), 



(27) 



with standard deviation a. By taking logs of equation[T4]we further 
see that 



15 



■ Simulation Values 

Prediction fi 4.250 
— Transition Region /3 = 4.2 50 ±0.250 




NO FRAGMENTATION - 
STABLE 




3.5 



Figure 5. Plot of /?„ at fragmentation for various values of n. The error 
bars correspond to the greatest non-fragmenting and smallest fragmenting 
values of fi found in the simulations, and the cross-hatched transition region 
represents uncertainty in the exact value of /?o- Note also that where n < 
discs may become thermally unstable. 



In Ofooi = In/? - n In T + n In T . 



(28) 



A standard property of the Normal distribution is that for a Nor- 
mally distributed random variable X ~ N(ji, a 2 ), the distribution of 
aX + bis given by N(afi + b, a 1 a 1 ). Hence from equation|28]we see 
that the distribution of In Q.t caol at fragmentation is such that 



In Qt c , 



A'GnjS^nV), 



(29) 



i.e., the distribution of In Q.t coo ] is centred around In/?,, for all n, 
reducing to a f5-function in the limit where n becomes zero and be- 
coming more spread out as n becomes large. Thus in order to coun- 
teract the increased width of the distribution, and thus the increased 
fraction of the gas that is below the fragmentation threshold, the av- 
erage must rise. This is clearly illustrated in Fig. [6] for values of n 
between and 4, and where f}„ is given in each case by equation[26l 
with/?o = 4.25. 



5 OPACITY-BASED ANALYTIC DISC MODELS 

Having quantified the effects of a temperature-dependent cooling 
law on the fragmentation boundary of protoplanetary discs, we are 
now in a position to use the known cooling laws for each opacity 
regime (as given by equationfTJt to determine the dominant cooling 
mechanisms throughout the radial range. We can therefore also use 
this to re-evaluate the regions of such discs that are unstable to 
fragmentation , in a similar manner to the analysis undertaken by 
IClarke (2009T) . 

In order to do this in a physically realistic manner we must 
also take into account the effects of the magneto-rotational in- 
stability (MRI), which operates when the disc becomes suffi- 
ciently ionised. Considering only thermal ionisation, we assume 
that the MRI becomes active when the disc temperature rises above 
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Figure 6. Variation in the distribution of In £2f coo l as a function of n, clearly 
showing the increasing width of the distribution with increasing n. Note that 
in the case where n = the distribution reduces to a iS-function. 



1000 K dClarke. 20091) . Although est imates of the vis cosity pro- 
vided through this instability vary (see lKing et al. (20071) for a sum- 
mary), numerical simul ations suggest it should be in the range 
0.001 < q-mri < 0.01 dWinters et al.. 200 jSano et al.. 2004 . We 
therefore assume that the MRI is the dominant instability in the 
disc wherever T > 1000K and the a delivered by the gravitational 
instability falls below 0.01. 

To obtain the disc temperature, we note that equations [311711 101 
QT]and[T2]self-consistently allow the disc properties to be evaluated 
for any given stellar mass M„ mass accretion rate M and radius R, 
when combined with the relation 

3act 

(30) 



GQ 



(see for instance lciarke 2009kafikov 2009bice & Armitage 20091) . 
We can thus derive the dependence of the disc temperature T on Q, 
M m , R and M, such that 



9k \GYRj \2nl 



(3D 



to prevent the temperature becom- 
we assume a fiducial background tem- 



Finally, in order 
ing too low, 

perature for the interstellar medium (ISM) of 10A" 
( ID' Alessio et al.. 199dHartmann et al., 1998h . In this case, we 
no longer assume that equation [3] holds, as there is additional 
heating from the background as well as from the gravitational 
instability. 

Since there is a strong dependence on temperature in certain 
opacity regimes (see Table [T} it is important that the equation of 
state adequately captures the correct behaviour of both the ratio 
of specific heats y and the mean molecular weight /j, as varia- 
tion in these can have significant effects on the system overall. To 
implement the equation of state we therefore make the assump- 
tion that the gas phase of the disc contains only hydrogen and he- 
lium, in the ratio 70: 30. We can make this assumption because al- 



Figure 7. Plot showing the variation of the ratio of specific heats y as a 
function of temperature and density. The density corresponds to the quoted 
radii for a Q = 1 disc about a 1 M Q star. Note that the inverse function for 
temperature in terms of y, T(y, p) is multi- valued. 



though the metallicity of the disc is important for the opacity (and 
thus the cooling), it makes very little contribution to the equation 
of state. Furthermore, the ratio of ortho- to para-hydrogen is as- 
sumed to be held constant at 3 : 1. Following on from the analysis of 
iBlack & Bodenheimer (1975t) . IStamatellos et al. (2007bl) produced 
tabulated values of p, T, y and fj for this equation of state and it 
is these values that we have used throughout. The variation of y 
with both temperature and density is shown in Fig. [7j - for the 
variation of the mean mo lecu lar weight p the rea der is referred to 
Stamatellos et al. (2007bl) and lForgan et al. (20091) . 

With this tabulated equation of state we can now solve the sys- 
tem of equations for Qt coo \ for any given values of Q, R, M and M, 
for each opacity regime. For simplicity we assume that the system 
is marginally gravitationally stable throughout, such that Q = I. 
Furthermore, since we know the dependence of Clt coo \ on tempera- 
ture for each of the opacity regimes, we can use equation |26l (with 
jSo = 4.25) to predict the (average) value of Q.t coo \ at which we would 
expect fragmentation, the results of which are shown in Table [4] 
Note that since the value of /S„ depends only on the relative size of 
the perturbations in temperature and not on either the mean temper- 
ature itself or the value of Q, we do not expect to see any variation 
in yS„ with varying Q, whereas the value of fl? coo | will vary with 
both. Using equations [3J] and [T3] we find that the Q dependence of 
flf coo i is 



-l+3(o+l)/(2i-7) 



(32) 



and thus except where b » 3.5 (such as in the regime where molec- 
ular line cooling dominates the opacity) the effects of Q variation 
are small. Nonetheless, in all optically thick cases, the effect of an 
increase in Q is to decrease the value of 0,t coo \, as can be seen from 
Table |3 

In Fig. [8] we therefore show the variation in Df coo i for a disc 
about a 1M G protostar as a function of radius at mass accretion 
rates of 10 -4 , 10" 6 , 10" 7 and 10" 8 M s yr~'. (For completeness, the 
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Figure 8. Value of Clt coo \ as a function of radius for accretion rates of 1CT 4 (top left), 1(T 6 (top right), 1(T 7 (bottom left) and 1CT 8 (bottom right) M Q yr , 
for a disc about a 1M Q star. The unshaded regions are optically thick (r > 5), the horizontally shaded areas are transitional (0.2 < r < 5) and the cross-hatched 
regions are optically thin (r < 0.2). The vertically shaded areas denote regions of the disc that are MRI active. The disc is stable against fragmentation 
wherever the value of flf CO ol is greater than the fragmentation boundary (shown by the heavy solid line). The dotted lines show the values that fi( coc j and the 
fragmentation boundary would take if the MRI were not active. 



various opacity regimes are shown in Fig.[9]for an accretion rate of 
10~ 4 M Q yr~' - all other accretion rates are qualitatively similar.) 
From the lower two panels (where the accretion rates are 1(T 7 and 
1CT 8 M Q yr~' for the left and right panels respectively) we see that 
at low accretion rates the fragmentation boundary becomes fixed at 
approximately 50AU, and that this is unaffected by the transition to 
the optically thin regime. This is down to the fact that the temper- 
ature becomes limited below by the background ISM temperature 
of 10K, and is therefore decoupled from the mass accretion rate. 

As the accretion rate rises to ~ 1CT 4 M G yr~' however, the 
disc becomes unstable to fragmentation at a wide range of radii 
due to the increase in the fragmentation boundary caused by the 
temperature dependence. Although an island of stability exists be- 
tween approximately 10 - 25 AU (where cooling is dominated by 
dust grains), all other radii become unstable. 

Note also that at low radii the disc becomes MRI active. This 
occurs at radii from ~ 1 - SAU dependent on M, which corresponds 
roughly to the transition to the dust sublimation opacity regime. 
For accretion rates of M < 10 4 M yr~' Fig. [8] suggests that the 
disc will be stable against fragmentation when the MRI is active, 
as in the absence of the MRI the value of fl/ coo i would be above the 
fragmentation boundary. However, where M a 1CT 4 M Q yr -1 the 



picture is less clear, as the disc is MRI active whilst simultaneously 
being unstable to fragmentation. However, iFromang et al. (20041) 
have suggested that where both instabilities operate the interaction 
causes the gravitationally-induced stress to weaken by a factor of 
two or so, which may stabilise the region against fragmentation. 

Nonetheless, throughout the range of mass accretion rates in- 
vestigated here there are no purely self-gravitating solutions at low 
radii, as the MRI is always active. It is however clear that for radii 
of ~ 5 - 50AU the susceptibility to fragmentation of a disc depends 
strongly on its steady state accretion rate, and that beyond approxi- 
mately 50AU, with a 10X" background temperature discs are always 
unstable to fragmentation. 

Finally it is useful to see how the fragmentation and MRI 
boundaries vary as a function of both R and M, and this is shown 
in Fig. [10] assuming that as before the central protostar has mass 
M, = 1M G . Here we have also included the fragmentation bound- 
ary in the case where the effects of temperature perturbations are ig- 
nored, i.e. where /? = 4.25 at fragmentation fo r all opacity re gimes, 
which allows for comparison with the work of lClarke (20091) . 

Fig. [10] shows clearly that by including the effects of temper- 
ature perturbations, the mass accretion rate at which fragmentation 
occurs is reduced, with an increased effect as the dependence of 
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Table 4. Predictions for the fragmentation boundary f}„ for each opacity 
regime in the optically thick case. The italicised case gives the prediction 
in the optically thin limit for ices, the only regime in our models where the 
disc becomes optically thin. Note that for large positive exponents (such 
as for hydrogen scattering) the value of /?„ becomes undefined. Note also 
that where the temperature exponent n is positive the regime may become 
susceptible to thermally instabilities. 



Opacity Regime 



Dependence of 0.t coa \ 
on T on Q 



Ices 


- none - 


or 1 


4.250 


Ices' 
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Q2/3 
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Figure 9. Plot of fiZ coo i for a mass accretion rate of 1 4 M Q yr 1 indicating 
the effects of the various opacity regimes. 



nf coo | on temperature increases. As before we note that there is now 
a region with M « 1CT 4 M G yr~' and R < 10AU where both the 
MRI is active and the disc is unstable to fragmentation. For accre- 
tion rates of ~ 10~ 5 - 1CT 3 M Q yr -1 there are limited radial ranges 
where a marginally gravitationally stable state exists, with regions 
that are unstable to fragmentation at both higher and lower radii. 

Fig. [TO] also shows how the stability of the disc to fragmenta- 
tion varies with the background ISM temperature. For low mass ac- 
cretion rates we see that as the background temperature decreases, 
the disc actually becomes stable out to larger radii. This can be ex- 
plained as follows: In the optically thin case where the cooling is 
dominated by ices (the regime in which this phenomenon is found) 
the value of Q.t CBO ] is given by 



fi/ c , 



1 



8<xk ji(y- 1) 



Q U2 p V2 T - 



(33) 




Figure 10. Plot showing the regions expected to be marginally gravitation- 
ally stable (unshaded), unstable to fragmentation (horizontal shading) and 
unstable to the magneto-rotational instability (vertical shading) in a disc 
about a 1M Q protostar. The cross-hatched regions show where where the 
disc is unstable to both the MRI and fragmentation. The more widely spaced 
horizontally shaded region to the lower right would become unstable to 
fragmentation if the minimum temperature limit of 10K was removed, and 
the fragmentation boundary moves to the right as the minimum tempera- 
ture is decreased. The short dashed line corres ponds to the f ragmentation 
boundary if a fixed value of/3 = 4.25 is used (cf. lClarke 2009b . 



37? VGMr 



1 



-R- 



(34) 



SCTK fj{j - 1) 

where we have used equation \\0\ to eliminate p in equation [34] 
Hence at a fixed radius R = i?f rag , increasing the temperature T 
decreases Q.t cao i and thereby destabilises the disc. Eventually, for 
some T = Tf mg we reach Qf coo i = 15.570 (from Tabled and the 
disc becomes unstable to fragmentation. 

From equation [34] we see that on the fragmentation boundary 
(where by construction Q.t coo i = 15.570 is constant), rf rae ~ R h ll'°- 
Now assuming that the temperature at which fragmentation occurs 
is at or above the background temperature (i.e. T^ag > r min ) then 
equation [30] holds, and we find similarly that the accretion rate at 
fragmentation Mf llle is given by Mf rag ~ . We therefore find that 
the radius as which fragmentation occurs increases with decreas- 
ing accretion rate such that Rf mg ~ M h ^ 19 . Hence, decreasing the 
background temperature decreases the accretion rate at which the 
disc becomes unstable to fragmentation, and likewise increases the 
radius at which this occurs. 

Note however that once T {mg is below the background tem- 
perature, (i.e. when T fmg < T min ) the disc temperature becomes de- 
coupled from the accretion rate, and hence all accretion rates below 
^min = ^fraa(^min) 316 unstable to fragmentation for radii R > /? frag . 



6 DISCUSSION AND CONCLUSIONS 

In summary, we have found from controlled numerical experiments 
with an imposed temperature dependent cooling law that the ef- 
fect of temperature dependence is to increase the value of fl/ coo i 
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at which the disc will fragment into bound objects. Furthermore, 
this tendency to fragment is greater the more strongly the cooling 
function depends on the local disc temperature. In t his respect, this 
confirms the results of Ijohnson & Gammie (20031) . who likewise 
noted a markedly increased tendency towards fragmentation in cer- 
tain opacity regimes. This result has been a ttributed to un certainty 
in the value of Q in the self-regulated state dClarke. 2009b . equiva- 
lent to uncertainty in the equilibrium temperature in our models. 

However, our results show that this is only one of two mecha- 
nisms that affect the fragmentation boundary, and one that we have 
been able to account for a posteriori by using effective values of f} 
rather than those input to the simulations. The other effect is due 
to the strength of the intrinsic temperature perturbations about the 
mean. In the case where the cooling law is dependent on temper- 
ature, perturbations about the equilibrium temperature will mean 
that some fraction of the gas has a lower value of f2f coo i than aver- 
age. Once this fraction reaches a critical value, the disc will become 
unstable to fragmentation. As the dependence of the cooling on 
these temperature perturbations increases, at a given average value 
of flf C00 ] the percentage of gas that lies below the critical value also 
increases, and thus the average must increase to avoid fragmenta- 
tion. 

We therefore find that the effect of allowing the cooling func- 
tion to depend on the local temperature is to make the disc more 
unstable to fragmentation, and we have been able to quantify this 
variation (see equation I26t . Combining this with predictions of 
the temperature dependence of protoplanetary discs using opacity- 
based cooling functions, we find that the fragmentation boundary 
can be increased by approximatel y an order of magnitude in t erms 
of Qfeooi, m close agreement with ljohnson & Gammie (2003b . We 
have also found that the RMS strength of the temperature perturba- 
tions can be correlated to the average cooling strength (see equa- 
tion|25ll, in a ver y similar manner to t hat found for the surface den- 
sity fluctuations dCossins et al.. 2009b . 

Using these predicted values in analytic models of marginally- 
gravitationally stable Q = 1 discs with a representative equation of 
state, we have found that the susceptibility of such discs to frag- 
mentation into bound objects is also sensitive to the steady state 
mass accretion rate, as shown in Fig. [TO] Others have noted that 
in the optically thick limit where the opacity is dominated by ices, 
flkooi is independent of temperature, and thus the cooling rate is 
determined only by the local density, itself a func tion of radius 
dMatzner & Levin. 2005feafikov. 200jClarke. 2009b . It has there- 
fore been suggested that once the cooling becomes dominated by 
ices fragmentation beyond some radius on the order of 100 A U be- 
comes inevitable, and indeed we find that with a background ISM 
temperature of 10A^, fragmentation occurs at ~ 50AU for all accre- 
tion rates below ~ 10~ 5 M Q yr -1 . 

However, if this minimum temperature condition is relaxed, 
we find that the change in cooling due to entering the optically 
thin regime has the effect of stabilising the disc out to large radii. 
(The fact that allowing it to become cooler actually stabilises the 
disc is due to the fact that in this regime Q.t coo \ increases with de- 
creasing temperature, and thus a hot disc has a shorter cooling 
time than a cold one.) For Class II / Classic T Tauri objects em- 
bedded in a cold medium with accretion rates below a few times 
10" 7 M G yr -1 , it is therefore possible that extended discs well be- 
yond 100 A U may be stable against fragmentation (they may well 
be stable against gravitational instabilitites altogether), and indeed 
discs with radii of at leas t 200 AU have been observed (see for ex- 
ample lEisner et ah 2008b . Nonetheless, discs with accretion rates a t 
the higher end of the scale (M ss 10~ 6 M G yr~ ' . lHartmann 2009b 



will still be unstable to fragmentation at radii beyond ~ 50AC7. It 
should be borne in mind however that in the outer regions of discs 
where the surface density is low, non-thermal ionisation (from cos- 
mic rays. X-rays etc) can trigger the MRI, and this may provide an 
alternative m echanism for preventing fragmentation, as shown in 
IClarke (20091) . 

Fig.[TOJalso shows another important result, that for accretion 
rates between 10~ 8 - 10" 2 M G yr~' discs cannot exist in a non- 
fragmenting purely self-gravitating state at radii < 2 - 5AU. In this 
regime discs are either MRI active (M < 10 -4 M Q yr ) or un- 
stable to fragmentation (M > 10~ 4 M e yr"" 1 ). We also find that 
in a narrow band of accretion rates ~ 10~ 4 M Q yr~' it is possible 
for discs to be both MRI active and unstable to fragmentation, al- 
thou gh the exact interact ion of these two instabilities is uncertain 
(see Fromang et al. 2004) . It is therefore the case that for steady- 
state protoplanetary discs the gravitational instability cannot drive 
accretion directly onto the protostar - either the MRI or the ther- 
mal instability m ust act at low radii, as has been prop osed for FU 
Orionis outbursts dArmitage et al. 200llZhu et al. 2009b . 

Finally, our results agree with the generally accepted 
view that planet formation through gravitationally-induced 
fragmentation is unlikely to occur at radii less than 50 - 100 AU 



:ely t 

dMatzner & Levin. 2005kafikov. 2005]Whitworth & Stamatellos. 200dClarke. 2009b 
although this critical radius varies with both the mass accretion 
rate and the background ISM temperature. Within this radius the 
core accretion model remains likely to be the dominant mode of 
planet formation. Outside this radius however, the fragmentation 
of spira l arms will pr oduce gaseous planets, a result which matches 
that of iBolev (2009b using a grid-based hydrodynamical model 
with radiative transfer - fragmentation was noted at ~ 100AI/ 
about a IM Q protostar. This r esult is further corroborated by 
IStamatellos & Whitworth (20 08) whose radiative transfer SPH 
code suggested a massive disc about a 0.7M o protostar would 
rapidly fragment into planetary mass objects or brown dwarf 
companions beyond approximately 100AI/. Although the mass 
accretion rate onto the central object is not stated in either case, we 
find that these figures are nonetheless in general agreement with 
our predictions. 
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